home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_10 / a10_2.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.7 KB  |  131 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 10.2 (Forward-Difference Method for the Heat Equation).
  9. % Section    10.2,    Parabolic Equations, Page 516
  10. echo on; clc; format short; hold off; clear
  11.  
  12. % PARABOLIC EQUATIONS.
  13. % Forward difference solution for the heat equation
  14.  
  15. %                          2
  16. %            u (x,t)   =  c  u  (x,t)
  17. %             t               xx
  18.  
  19. %  u(x,0) = f(x)   for  0 < x < a     and
  20.  
  21. %  u(0,t) = g1(t) and u(a,t) = g2(t)  for  0 ≤ t ≤ b
  22.  
  23. % A numerical approximation is computed over the rectangle
  24.  
  25. %         0 ≤ x ≤ a ,   0 ≤ t ≤ b.
  26.  
  27. pause    % Press any key to continue.
  28.  
  29. clc;
  30. %    Store f(x), g1(x), g2(x) in  f.m  g1.m  g2.m
  31. % function z = f(x)
  32. % z = 4*x - 4*x^2;
  33.  
  34. % function z = g1(t)
  35. % z = 0;
  36.  
  37. % function z = g2(t)
  38. % z = 0;
  39.  
  40. delete f.m
  41. diary f.m; disp('function z = f(x)');...
  42.            disp('z = 4*x - 4*x^2;');...
  43. diary off;
  44.  
  45. delete g1.m
  46. diary g1.m; disp('function z = g1(t)');...
  47.             disp('z = 0;');...
  48. diary off;
  49.  
  50. delete g2.m
  51. diary g2.m; disp('function z = g2(t)');...
  52.             disp('z = 0;');...
  53. diary off;
  54.  
  55. % Remark. f.m g1.m g2.m forwdif.m are used for Algorithm 10.2
  56. f(0); g1(0); g2(0); 
  57. pause    % Press any key to continue.
  58.  
  59. clc;
  60. %    Place the endpoint of [0,a] in  a.
  61.  
  62. %    Place the endpoint of [0,b] in  b.
  63.  
  64. %    For the heat equation,  enter the constant  c.
  65.  
  66. %    Over [0,a] enter the number of grid points  n.
  67.  
  68. %    Over [0,b] enter the number of grid points  m.
  69.  
  70. a  =  1.0;
  71. b  =  0.2;
  72. c  =  1;
  73. n  =  6;
  74. m  =  11;
  75.  
  76. U = forwdif('f','g1','g2',a,b,c,n,m);
  77.  
  78. pause    % Press any key to see the solution. 
  79.  
  80. clc; clg;
  81. Z = fliplr(U);
  82. mesh(Z);...
  83. Mx1 = 'The forward difference solution to the heat equation.';...
  84. title(Mx1);...
  85. shg; pause % Press any key to continue.
  86.  
  87. W = U';
  88. points = W(:,2:n-1);
  89. Mx2 = 'The forward difference method is stable because r <= 1/2.';
  90. clc,echo off,diary output,disp(' '),...
  91. disp(Mx1),disp(' '),disp(points),disp(Mx2),...
  92. diary off,echo on
  93.  
  94.  
  95. pause    % Press any key to continue.
  96.  
  97. clc;
  98. %    Place the endpoint of [0,a] in  a.
  99.  
  100. %    Place the endpoint of [0,b] in  b.
  101.  
  102. %    For the heat equation,  enter the constant  c.
  103.  
  104. %    Over [0,a] enter the number of grid points  n.
  105.  
  106. %    Over [0,b] enter the number of grid points  m.
  107.  
  108. a  =  1.0;
  109. b  =  0.333333333;
  110. c  =  1;
  111. n  =  6;
  112. m  =  11;
  113.  
  114. U = forwdif('f','g1','g2',a,b,c,n,m);
  115.  
  116. pause    % Press any key to see the solution. 
  117.  
  118. clc; clg;
  119. Z = fliplr(U);
  120. mesh(Z);
  121. title(Mx1)
  122.  
  123. W = U';
  124. points = W(:,2:n-1);
  125. Mx3 = 'The forward difference method is unstable because r > 1/2.';
  126. clc,echo off,diary output,disp(' '),...
  127. disp(Mx1),disp(' '),disp(points),disp(Mx3),...
  128. diary off,echo on
  129.  
  130.  
  131.